Objectives:
# open raster data
lidar_dem <- raster(x = "data/week-03/BLDR_LeeHill/pre-flood/lidar/pre_DTM.tif")
# plot raster data
plot(lidar_dem,
main = "Digital Elevation Model - Pre 2013 Flood")
# zoom in to one region of the raster
plot(lidar_dem,
xlim = c(473000, 473030), # define the x limits
ylim = c(4434000, 4434030), # define y limits for the plot
main = "Lidar Raster - Zoomed into one small region")
# view resolution units
crs(lidar_dem)
## CRS arguments:
## +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
# assign crs to an object (class) to use for reprojection and other tasks
myCRS <- crs(lidar_dem)
myCRS
## CRS arguments:
## +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
# what is the x and y resolution for your raster data?
xres(lidar_dem)
## [1] 1
yres(lidar_dem)
## [1] 1
# plot histogram
hist(lidar_dem,
main = "Distribution of surface elevation values",
# breaks = c(1600, 1800, 2000, 2100),
# breaks = 3,
xlab = "Elevation (meters)", ylab = "Frequency",
col = "purple")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 1% of the raster cells were used. 100000 values used.
lidar_dsm <- raster(x = "data/week-03/BLDR_LeeHill/pre-flood/lidar/pre_DSM_hill.tif")
# plot raster data
plot(lidar_dsm,
main = "Digital Surface Model")
# plot histogram
hist(lidar_dsm,
main = "Digital Surface Model - Histogram",
# breaks = c(1600, 1800, 2000, 2100),
# breaks = 3,
xlab = "Elevation", ylab = "Frequency",
col = "purple")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 1% of the raster cells were used. 100000 values used.
A GeoTIFF is a standard .tif or image file format that includes additional spatial (georeferencing) information embedded in the .tif file as tags. These are called embedded tags, tif tags. These tags can include the following raster metadata:
# view attributes associated with your DTM geotiff
GDALinfo("data/week-03/BLDR_LeeHill/pre-flood/lidar/pre_DTM.tif")
## Warning in GDALinfo("data/week-03/BLDR_LeeHill/pre-flood/lidar/
## pre_DTM.tif"): statistics not supported by this driver
## rows 2000
## columns 4000
## bands 1
## lower left origin.x 472000
## lower left origin.y 4434000
## res.x 1
## res.y 1
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs
## file data/week-03/BLDR_LeeHill/pre-flood/lidar/pre_DTM.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float32 TRUE -3.402823e+38 128 128
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 -4294967295 4294967295 NA NA
## Metadata:
## AREA_OR_POINT=Area
# view attributes / metadata of raster
# view crs
crs(lidar_dem)
## CRS arguments:
## +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
lidar_dem@extent
## class : Extent
## xmin : 472000
## xmax : 476000
## ymin : 4434000
## ymax : 4436000
extent_lidar_dsm <- extent(lidar_dsm)
extent_lidar_dem <- extent(lidar_dem)
# Do the two datasets cover the same spatial extents?
if (extent_lidar_dem == extent_lidar_dsm) {
print("Both datasets cover the same spatial extent")
}
## [1] "Both datasets cover the same spatial extent"
compareRaster(lidar_dsm, lidar_dem, extent = TRUE)
## [1] TRUE
compareRaster(lidar_dsm, lidar_dem, res = TRUE)
## [1] TRUE
# single layer (or band) vs multi-layer (Band Geotiffs)
nlayers(lidar_dsm)
## [1] 1
DSM - DTM = CHM
# plot raster data
plot(lidar_dem,
main = "Digital Elevation Model - Pre 2013 Flood")
# open raster data
lidar_dsm <- raster(x = "data/week-03/BLDR_LeeHill/pre-flood/lidar/pre_DSM.tif")
# plot raster data
plot(lidar_dsm,
main = "Lidar Digital Surface Model (DSM)")
# open raster data
lidar_chm <- lidar_dsm - lidar_dem
# plot raster data
plot(lidar_chm,
main = "Lidar Canopy Height Model (CHM)")
# plot raster data
plot(lidar_chm,
breaks = c(0, 2, 10, 20, 30),
main = "Lidar Canopy Height Model",
col = c("white", "brown", "springgreen", "darkgreen"))
# check to see if an output directory exists
dir.exists("data/week-03/outputs")
## [1] TRUE
# if the output directory doesn't exist, create it
if (!dir.exists("data/week-03/outputs")) {
# if the directory doesn't exist, create it
# recursive tells R to create the entire directory path (data/week-03/outputs)
dir.create("data/week-03/outputs", recursive=TRUE)
}
# # export CHM object to new GeotIFF
# writeRaster(lidar_chm, "data/week-03/outputs/lidar_chm.tiff",
# format = "GTiff", # output format = GeoTIFF
# overwrite = TRUE) # CAUTION: if this is true, it will overwrite an existing file
# Challenge: Calculate changes in terrain
lidar_dem_post <- raster(x = "data/week-03/BLDR_LeeHill/post-flood/lidar/post_DTM.tif")
lidar_dem_diff <- lidar_dem_post - lidar_dem
plot(lidar_dem_diff,
main = "DEM difference in elevation before and after flood")
lidar_dsm_post <- raster(x = "data/week-03/BLDR_LeeHill/post-flood/lidar/post_DSM.tif")
lidar_chm_post <- lidar_dsm_post - lidar_dem_post
lidar_chm_diff <- lidar_chm_post - lidar_chm
plot(lidar_chm_diff,
main = "CHM difference in elevation before and after flood")
# writeRaster(lidar_chm_diff, "data/week-03/outputs/lidar_chm_diff.tiff",
# format = "GTiff", # output format = GeoTIFF
# overwrite = TRUE) # CAUTION: if this is true, it will overwrite an existing file
# open canopy height model
lidar_chm <- raster("data/week-03/BLDR_LeeHill/outputs/lidar_chm.tif")
summary(lidar_chm)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (1.25% of all cells)
## lidar_chm
## Min. 0.0000000
## 1st Qu. 0.0000000
## Median 0.0000000
## 3rd Qu. 0.6900635
## Max. 23.8999023
## NA's 0.0000000
# plot histogram of data
hist(lidar_chm,
xlim = c(2, 25),
ylim = c(0, 4000),
main = "Distribution of raster cell values in the DTM difference data",
xlab = "Height (m)", ylab = "Number of Pixels",
col = "springgreen")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 1% of the raster cells were used. 100000 values used.
# see how R is breaking up the data
histinfo <- hist(lidar_chm)
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 1% of the raster cells were used. 100000 values used.
histinfo$counts
## [1] 76109 3462 3041 2880 2468 2161 1925 1777 1512 1229 927
## [12] 711 588 436 293 190 114 82 42 28 16 2
## [23] 5 2
histinfo$breaks
## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
## [24] 23 24
# zoom in on x and y axis
hist(lidar_chm,
xlim = c(2, 25),
ylim = c(0, 1000),
breaks = 100,
main = "Histogram of canopy height model differences \nZoomed in to -2 to 2 on the x axis",
col = "springgreen",
xlab = "Pixel value")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 1% of the raster cells were used. 100000 values used.
# You may want to explore breaks in your histogram before plotting your data
# Note: When using breaks = c() in the plotting, the y axis will auto change to density
# density = (frequency/total samples)/bin width
lidar_info <- hist(lidar_chm,
breaks = c(0, 2, 4, 7, 30),
ylim = c(0,1),
main = "Histogram with custom breaks",
xlab = "Height (m)" , ylab = "density",
col = "springgreen")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 1% of the raster cells were used. 100000 values used.
lidar_info$breaks
## [1] 0 2 4 7 30
lidar_info$counts
## [1] 79484 5982 6615 7919
lidar_info$density
## [1] 0.397420000 0.029910000 0.022050000 0.003443043
# create classification matrix
reclass_df <- c(0, 2, NA,
2, 4, 1,
4, 7, 2,
7, Inf, 3)
reclass_df
## [1] 0 2 NA 2 4 1 4 7 2 7 Inf 3
# reshape the object into a matrix with columns and rows
reclass_m <- matrix(reclass_df,
ncol = 3,
byrow = TRUE)
reclass_m
## [,1] [,2] [,3]
## [1,] 0 2 NA
## [2,] 2 4 1
## [3,] 4 7 2
## [4,] 7 Inf 3
# reclassify the raster using the reclass object - reclass_m
chm_classified <- reclassify(lidar_chm,
reclass_m)
# view reclassified data
barplot(chm_classified,
main = "Number of pixels in each class")
## Warning in .local(height, ...): a sample of 12.5% of the raster cells were
## used to estimate frequencies
# str(chm_classified)
# assign all pixels that equal 0 to NA or no data value
chm_classified[chm_classified == 0] <- NA
# create color object with nice new colors!
chm_colors <- c("palegoldenrod", "palegreen2", "palegreen4")
# plot reclassified data
plot(chm_classified,
legend = FALSE,
col = chm_colors,
axes = FALSE, # remove the box around the plot
box = FALSE,
main = "Classified Canopy Height Model \n short, medium, tall trees")
legend("topright",
legend = c("short trees", "medium trees", "tall trees"),
fill = chm_colors,
border = FALSE,
bty = "n")
# import the vector boundary
crop_extent <- readOGR("data/week-03/BLDR_LeeHill/clip-extent.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "data/week-03/BLDR_LeeHill/clip-extent.shp", layer: "clip-extent"
## with 1 features
## It has 1 fields
## Integer64 fields read as strings: id
# plot imported shapefile
# notice that you use add = T to add a layer on top of an existing plot in R.
plot(crop_extent,
main = "Shapefile imported into R - crop extent",
axes = TRUE,
border = "blue")
# crop the lidar raster using the vector extent
lidar_chm_crop <- crop(lidar_chm, crop_extent)
plot(lidar_chm_crop, main = "Cropped lidar chm")
# add shapefile on top of the existing raster
plot(crop_extent, add = TRUE)
# crop_extent <- readOGR("data/week-03/BLDR_LeeHill/clip-extent.shp")
lidar_chm_diff_crop <- crop(lidar_chm_diff, crop_extent)
# create classification matrix
reclass_df <- c(-25, 0, 1,
0, 25, 2)
reclass_df
## [1] -25 0 1 0 25 2
# reshape the object into a matrix with columns and rows
reclass_m <- matrix(reclass_df,
ncol = 3,
byrow = TRUE)
reclass_m
## [,1] [,2] [,3]
## [1,] -25 0 1
## [2,] 0 25 2
# reclassify the raster using the reclass object - reclass_m
chm_classified <- reclassify(lidar_chm_diff_crop,
reclass_m)
# view reclassified data
barplot(chm_classified,
main = "Number of pixels in each class")
## Warning in .local(height, ...): a sample of 14.3% of the raster cells were
## used to estimate frequencies
# str(chm_classified)
# assign all pixels that equal 0 to NA or no data value
chm_classified[chm_classified == 0] <- NA
# create color object with nice new colors!
chm_colors <- c("blueviolet", "aquamarine2")
# plot reclassified data
plot(chm_classified,
legend = FALSE,
col = chm_colors,
axes = FALSE, # remove the box around the plot
box = FALSE,
main = "Classified Canopy Height Change Model")
legend("bottomright",
legend = c("negative change", "positive change"),
fill = chm_colors,
border = FALSE,
bty = "n")
# crop_extent <- readOGR("data/week-03/BLDR_LeeHill/clip-extent.shp")
lidar_dem_diff_crop <- crop(lidar_dem_diff, crop_extent)
# create classification matrix
reclass_df <- c(-20, 0, 1,
0, 48, 2)
reclass_df
## [1] -20 0 1 0 48 2
# reshape the object into a matrix with columns and rows
reclass_m <- matrix(reclass_df,
ncol = 3,
byrow = TRUE)
reclass_m
## [,1] [,2] [,3]
## [1,] -20 0 1
## [2,] 0 48 2
# reclassify the raster using the reclass object - reclass_m
chm_classified <- reclassify(lidar_dem_diff_crop,
reclass_m)
# view reclassified data
barplot(chm_classified,
main = "Number of pixels in each class")
## Warning in .local(height, ...): a sample of 14.3% of the raster cells were
## used to estimate frequencies
# str(chm_classified)
# assign all pixels that equal 0 to NA or no data value
chm_classified[chm_classified == 0] <- NA
# create color object with nice new colors!
chm_colors <- c("blueviolet", "aquamarine2")
# plot reclassified data
plot(chm_classified,
legend = FALSE,
col = chm_colors,
axes = FALSE, # remove the box around the plot
box = FALSE,
main = "Classified Terrain Change Model")
legend("topright",
legend = c("negative change", "positive change"),
fill = chm_colors,
border = FALSE,
bty = "n")
myMap <- get_map(location = "Boulder, Colorado",
source = "google",
maptype = "terrain", crop = FALSE,
zoom = 6)
# plot map
ggmap(myMap)
myMap <- get_map(location = "Boulder, Colorado",
source = "google",
maptype = "satellite", crop = FALSE,
zoom = 6)
# plot map
ggmap(myMap)
# water color maps
myMap <- get_map(location = "Boulder, Colorado",
source = "stamen",
maptype = "watercolor", crop = FALSE,
zoom = 4)
# plot map
ggmap(myMap)
# add points to your map
# creating a sample data.frame with your lat/lon points
gage_location <- data.frame(lon = c(-105.178333), lat = c(40.051667))
# create a map with a point location for boulder.
ggmap(myMap) + labs(x = "", y = "") +
geom_point(data = gage_location, aes(x = lon, y = lat, fill = "red", alpha = 0.2), size = 5, shape = 19) +
guides(fill = FALSE, alpha = FALSE, size = FALSE)
# package 'maps'
map('state')
# add a title to your map
title('Map of the United States')
map('state', col = "darkgray",
fill = TRUE,
border = "white")
# add a title to your map
title('Map of the United States')
map('county', regions = "Colorado", col = "darkgray", fill = TRUE, border = "grey80")
map('state', regions = "Colorado", col = "black", add = TRUE)
# add the x, y location of the stream guage using the points
# notice i used two colors adn sized to may the symbol look a little brighter
points(x = -105.178333, y = 40.051667, pch = 21, col = "violetred4", cex = 2)
points(x = -105.178333, y = 40.051667, pch = 8, col = "white", cex = 1.3)
# add a title to your map
title('County Map of Colorado\nStream gage location')
map('state', fill = TRUE, col = "darkgray", border = "white", lwd = 1)
map(database = "usa", lwd = 1, add = TRUE)
# add the adjacent parts of the US; can't forget my homeland
map("state", "colorado", col = "springgreen",
lwd = 1, fill = TRUE, add = TRUE)
# add gage location
title("Stream gage location\nBoulder, Colorado")
# add the x, y location of hte stream guage using the points
points(x = -105.178333, y = 40.051667, pch = 8, col = "red", cex = 1.3)
# open dem hillshade
lidar_dem_hill <- raster(x = "data/week-03/BLDR_LeeHill/pre-flood/lidar/pre_DTM_hill.tif")
# plot raster data
plot(lidar_dem_hill,
main = "Lidar Digital Elevation Model (DEM)\n overlayed on top of a hillshade",
col = grey(1:100/100),
legend = FALSE)
plot(lidar_dem,
main = "Lidar Digital Elevation Model (DEM)",
add = TRUE, alpha = .5)